% This function computes the standard errors for the GMM estimates in
% the VAR(1) model when accounting for measurement errors in the factors

function out = getAVarTheta2(theta2,VarFactors,Cov10Factors,Cov01Factors,factors,bandwidthT)

% Computing the Jacobian matrix
epsValue            = 1e-9;
namesTheta2         = fieldnames(theta2);
numTheta2           = size(namesTheta2,1);
moments             = Var1_GMM_moments(theta2,VarFactors,Cov10Factors,Cov01Factors,factors);
T                   = size(moments,2);
%Q = mean(moments,2)'*mean(moments,2)
R                   = zeros(numTheta2,numTheta2);
for i=1:numTheta2
    theta2Grad = theta2;
    theta2Grad.(namesTheta2{i}) = theta2Grad.(namesTheta2{i}) + epsValue;
    moments_p_eps = Var1_GMM_moments(theta2Grad,VarFactors,Cov10Factors,Cov01Factors,factors);

    theta2Grad = theta2;
    theta2Grad.(namesTheta2{i}) = theta2Grad.(namesTheta2{i}) - epsValue;
    moments_m_eps = Var1_GMM_moments(theta2Grad,VarFactors,Cov10Factors,Cov01Factors,factors);
    
    R(i,:) = (mean(moments_p_eps(:,1:T-1),2)-mean(moments_m_eps(:,1:T-1),2))'/(2*epsValue);
end

% The Newey-West estimate of S with lags controlled by bandwidthT (Robust to auto-correlation and
% heteroskedasticity)
GAMMA0 = zeros(numTheta2,numTheta2);
for t=1:T-1
    GAMMA0 = GAMMA0 + moments(:,t)*moments(:,t)'/(T-1);
end
S      = GAMMA0;
GAMMAk = zeros(numTheta2,numTheta2,bandwidthT);
for k=1:bandwidthT
    for t=1+k:T-1
        GAMMAk(:,:,k) = GAMMAk(:,:,k) + moments(:,t)*moments(:,t-k)';
    end    
    GAMMAk(:,:,k) = GAMMAk(:,:,k)/(T-1-k);
    S = S + (1-k/(1+bandwidthT))*(GAMMAk(:,:,k) + GAMMAk(:,:,k)');
end

% The asympotic covariance matrix 
out.VarRobust  = 1/(T-1)*inv(R*inv(S)*R');
out.names      = namesTheta2;
